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Kinetic equations for the stationary state distribution function of ions moving through narrow 
pores are solved for a number of one-dimensional models of single ion transport. Ions move through 
pores of length L, under the action of a constant external field and of a concentration gradient. The 
interaction of single ions with the confining pore surface and with water molecules inside the pore 
are modelled by a Fokker-Planck term in the kinetic equation, or by uncorrelated collisions with 
thermalizing centres distributed along the pore. The temporary binding of ions to polar residues 
lining the pore is modelled by stopping traps or energy barriers. Analytic expressions for the 
stationary ion current through the pore are derived for several versions of the model, as functions 
of key physical parameters. In all cases, saturation of the current at high fields is predicted. Such 
simple models, for which results are analytic, may prove useful in the study of the current/voltage 
relations of ion channels through membranes. 

PACS numbers: 05.20.Dd, 51.10.+y, 87.16. Ac 



I. INTRODUCTION 



The flow of fluids through porous media is a classic 
problem which has many scientific and industrial applica- 
tions. For very narrow pores, with diameters of the order 
of 1 nm or less, continuum descriptions become inappli- 
cable and the transport of matter must be examined on 
the molecular scale. Examples include molecular or ionic 
permeation of zeolites [l| , of carbon nanotubes Q, y, Q 
and of aquaporins [j| and ion channels [|| through cell 
membranes. The simple kinetic models examined in this 
work are meant to represent ion channels: they are, how- 
ever, also more widely applicable - for example, we shall 
present results for ions flowing through an infinitely long 
pore that might represent a carbon nanotube or part of 
a zeolite. 

Ion channels arc pores through cell membranes, 
through which ions are transported under the influence 
of a concentration gradient and a large electric field. The 
permeability of the pores is highly selective for particular 
ions and the pores can also open and close to ion trans- 
port (a phenomenon known as "gating") in response to 
factors such as ligand binding or changes in the electric 
field or the membrane tension. Many channels contain 
a narrow region, the "selectivity filter" , where ionic mo- 
tion is essentially single-file 0, H, 0> LUl< Some chan- 
nels appear to transport only one ion at a time, while 
others use transport mechanisms involving multiple ions 
El U fill Il2| . Measurements of the current through 
individual channels have been possible for some time, 
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and these have resulted in a large amount of data, both 
on the gating characteristics of channels, and on their 
properties in the open state. These properties include 
the relationship between the ionic current and the elec- 
tric field applied across the membrane (current-voltage 
relations), as well as the conductivity of the channels 
as a function of the ionic concentration difference, at 
fixed applied voltage (conductance-concentration rela- 
tions). One of the challenges for theoreticians is to relate 
these functional characteristics to the geometric, phys- 
ical and chemical structure of the pores, which are be- 
coming increasingly well-known [l3|. This goal may be 
achieved by detailed simulations of the motion of ions 
and molecules th roug h s pec ific pores 0, Il5| . or using 
simplified models [la. Il7l Il8t Il9l ] . Rates of ion transport 
can be predicted directly, or by application of barrier- 
crossing theories such as Kramers rate theory [23, El • 
An alternative approach is the extension of continuum 
theories to the nanoscale. Goldman [23j and Hodgkin 
and Katz |24j](GHK), in their classic work, applied the 
one-dimensional diffusion equation in a constant electric 
field to predict current-voltage relations for ion channels. 
This work is generalized to specific and multi-dimensional 
ion channel models in the Poisson-Nernst-Planck (PNP) 
theory of ion channels [2^, [2(| , where numerical solution 
methods are used to obtain the current due to diffusion in 
the presence of complicated and self-consistent potential 
fields. 

In this paper, we explore the possibility of applying 
simple, analytically solvable, kinetic models to the prob- 
lem of transport in nanopores. Our approach is for the 
moment very general. We consider the motion of single 
ions through a one-dimensional pore of length L connect- 
ing two reservoirs at different ion concentrations, under 
the action of a constant electric field (we shall also con- 
sider the case of an infinitely long pore). Ions are ex- 
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pected to experience friction due to collisions with the 
inside surface of the pore and with other particles, such 
as water molecules, inside the pore. We include this ef- 
fect initially by a Langevin-like friction within the frame- 
work of a Fokker-Planck (FP) description, although we 
shall see that the FP approach presents some difficul- 
ties in the confined geometry of the finite length pore, 
which we shall attempt to overcome by using an alter- 
native description of the friction in terms of localized 
"thermalizing centres" . We also consider that there may 
be more specific binding interactions between the ion and 
the pore surface, for example with polar residues lining 
the surface. These arc modelled by "stopping traps" - on 
encountering such a trap, an ions is stopped, and later 
released to continue its motion under the influence of 
friction and the electric field. The action of the stop- 
ping traps on the ion may or may not depend upon its 
velocity. In all cases, we attempt to find general ana- 
lytic solutions for the stationary state ionic current j, and 
sometimes also for the ion distribution function f(x,v) 
(defined so that the probability of finding an ion between 
position x and x + dx with velocity between v and v + dv 
is f(x,v)dxdv). These solutions are functions of the ap- 
plied field, as well as of parameters such as the channel 
length, friction coefficient and probability density of stop- 
ping traps. We hope that these results may ultimately be 
used to analyse the transport behaviour of specific pores 
or channels, by adapting the above physical parameters 
to the known structure of the pore under consideration. 

The general kinetic equation for the model is presented 
in section ITU The kinetic equation is solved in section ITTTI 
for the case of a finite length channel with stopping traps 
but without friction, and in section IIVI in the case of a 
finite channel with friction but without traps; difficulties 
arising from the use of the FP operator in a pore of fi- 
nite length arc discussed. These difficulties do not arise 
in the case of an infinite (L — > oo) channel, for which a 
general solution is obtained in the presence of both fric- 
tion and traps, in section Returning to a finite L 
channel in section IVII the FP friction is replaced by a 
distribution of thermalizing centres throughout the pore, 
and an explicit expression for the current is obtained in 
the presence of such thermalizing centres and stopping 
traps. Concluding remarks are made in the last section. 



II. MODEL AND KINETIC EQUATION 

Our model channel is pictured in figure ^ The chan- 
nel, of length L, is located along the x-axis (-L/2 < 
x < L/2). The radius of the pore (which is assumed to 
be cylindrical) matches the ion radius, so that ionic mo- 
tion inside the pore is strictly one-dimensional. The pore 
links two reservoirs containing ionic solutions of linear 
concentrations pi (to the left) and p r (to the right): pi 
and p r are related to the bulk concentrations c/ and c r 
in the reservoirs by p r ^i = irR 2 c r j, where R is the radius 
of the pore. The inner surface of the pore is lined with 



Pt 



o 




a = q/E 

FIG. 1: Schematic view of the model channel. 



stopping traps of local average density p(x): we shall as- 
sume that the probability of finding n traps within the 
interval x\ < x < X2 is given by the Poisson distribution 
with parameter J x ^ p{x)dx. 

If an ion encounters a trap, its velocity is set to zero, 
generally irrespective of its initial velocity (although we 
shall also consider in section ITTTI the case of traps which 
discriminate between ions according to velocity). Inside 
the pore, the ion (of charge q and mass m) is subjected 
to a uniform electric field E, and hence undergoes an ac- 
celeration, towards the right, a = qE/m. After being 
stopped by a trap, the ion is therefore re-accelerated by 
the electric field. Ions also experience friction: this will 
initially be modelled by a force — jv, where v is the veloc- 
ity of the ion, as well as the thermalizing effect of a ran- 
dom force, although an alternative to this Langevin-like 
model will be presented in section IVTI In summary, the 
ion undergoes a constant acceleration due to the electric 
field, is slowed down by collisions with molecules inside 
the pore or on the pore surface (these processes being 
described by a friction process) and may be captured by 
traps along the channel, to account for temporary bind- 
ing to polar residues on the pore surface. We shall present 
results for stationary state ion flow only. 

The general kinetic equation for the stationary state 
ion distribution function f(x,v), in the presence of 
Poisson-distributed stopping traps of average density 
p(x) as well as a Langevin-like friction mechanism, with 
friction coefficient 7, is: 
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dw \w\f(x,w) - \v\f(x,v) 
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The left-hand side (l.h.s.) of (IJ describes free flow of 
ions under the action of the constant acceleration a aris- 
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ing from the external field. The r.h.s. contains two col- 
lision terms. The first accounts for the stopping traps: 
it is a balance between gain (in the population of zero 
velocity particles) and loss (of particles with velocity v). 
The second term is the Fokker-Planck operator acting 
on the distribution function: it accounts for the effect of 
the frictional and random forces. Note that the kinetic 
equation is for a single ion: it does not account for 
interactions between several ions within the pore. This 
limitation will be addressed in later work. 

The reservoirs on the left (x < —L/2) and on the right 
(x > L/2) of the channel are assumed to contain ions 
in thermodynamic equilibrium at the same temperature 
T, but generally at different densities: pi (to the left) 
and p r (to the right). The ion distribution functions in 
the reservoirs (not including the contribution of any ions 
coming out of the pore) are hence: 

fl(x,v) = pi4> T (v) fr(x,v) = Pr4* T (v) (2) 



where 



4> T (v) 



2nk B T 



exp 



2k B T 



(3) 



is the Maxwell velocity distribution function. 

For illustrative purposes, we first consider the case 
where acceleration, traps and friction are all absent, and 
an ion which enters the pore at one end keeps the same 
velocity until it reaches the other end. The ion distribu- 
tion function within the pore is then simply: 

f(x,v) = [p l 6(v)+p r 6(-v)]<j> T (v) (4) 



where 9 denotes the Hcaviside step function, 
current is given by: 



The 



j(x) = / dvvf(x,v) 



(5) 



For a stationary state, continuity requires that the cur- 
rent be independent of position: 



dj(x) 
dx 



= 



(6) 



Substituting J3J into JSJ one finds the result in the ab- 
sence of acceleration, traps or friction: 
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while the number density inside the channel 

/>oo 

n ( x ) — I dv f(x, v) 



(7) 



(8) 



is in this case given by n = (pi +p r )/2. Note that discon- 
tinuities arise in n(x) at the pore boundaries, x = ±L/2: 
this reflects the fact that the regions close to the pore 
mouth are not modelled in detail in this simple theory. 

In the subsequent sections, analytic solutions of the 
kinetic equation will be derived for the limiting cases 
7 = ( section llllll . p(x) = fsection llV|l and L — > oo 
(section 0). 



III. FINITE CHANNEL WITH TRAPS 

Consider a pore of finite length L, containing stopping 
traps but no friction mechanism. The traps have aver- 
age local density p(x) and are distributed according to a 
Poisson law as described in Section [H] on encountering 
such a trap, the velocity of an ion is reduced to zero. In 
the absence of friction, the kinetic equation for the 
stationary state ion distribution function f(x,v) simpli- 
fies to: 



(9) 



p(x) { S(v) 



-Hoc- 



dw\w\f{x,w) - \v\f(x,v) 



Note that we obtain Equation © (constant current 
throughout the pore) on integrating both sides of 10 over 
all velocities — oo < v < +oo. If the traps are on aver- 
age uniformly distributed (p(x) = p), Equation Q can 
be solved exactly for f(x,v) as shown in Appendix IaI 
However, an expression for the ionic current j can be ob- 
tained for any p(x) using simple arguments, without the 
need for an explicit solution for f(x,v). 

We first note that in the stationary state, the contribu- 
tion to the current due to an ion which enters the channel 
at one end depends only on its incoming velocity and its 
probability of eventually arriving at the other end, since 
j does not depend on x (and there are no interactions be- 
tween ions) . For this model, all ions entering the channel 
from the left reservoir at x — —L/2 will eventually reach 
x = L/2, since on being stopped by a trap they are re- 
accelerated by the field towards the right (assuming a 
is positive). Thus the contribution of these ions to the 
current is: 



]i = Pi 



k B T 
2ixm 



(10) 



Ions entering the channel from the right at x — L/2 will 
reach x = —L/2 only if they are not stopped either by 
the opposing field or by an encouter with a trap. An ion 
which is stopped is re-accelerated towards the right, so 
that it will exit the channel at L/2. The Poisson prob- 
ability of encountering no traps between X\ and x-i is: 



P(xi,x 2 ) = exp 



p(x)dx 



X\ < x 2 (11) 



In order to overcome the opposing field, ions must enter 
the channel with kinetic energy mv 2 /2 > maL, so that, 
assuming a Maxwell distribution of velocities at L/2, the 
distribution function f r (x,v) of particles which entered 
the channel at L/2 and which will eventually reach —L/2 
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f r (x,v)= p r 0(- V )9[- + 

exp 



2Trk B T 



L 

2 ' 
rn 



aL 



(12) 



k B T 



, L L 

xP , - 

' 2' 2 



from which Equation @ leads to a contribution to the 
current: 



Jr 



-Pr 



k B T 



maL 



exp 



L/2 



2nm ~" r ^ fcfiT J-l/2 
Adding I|1U|) and (|13|l leads to the result for j 

maL r L / 2 



p(x')dx' } (13) 



J = 



k B T 
2irm 



pi - p r exp 



k B T 



p(x')dx' 



-L/2 



For a uniform distribution of traps, Equation (|14|) re- 
duces to: 



J 



k B T 
27rm 



pi - p r exp < —L p 



k B T 



(15) 



in agreement with the full solution derived in Appendix 
lAl Equation l|14fl shows that in this model, in the absence 
of friction, the current does not depend on the spatial dis- 
tribution of the stopping traps, but only on the integral 
of p(x) between x = —L/2 and x = L/2. Note that 
Equations (|14fl and l|15fl were derived for a > 0. When 
a < 0, the roles of the right and left reservoirs must be 
interchanged. 

If the reservoir densities pi and p r are measured rel- 
ative to an arbitrary "reference density" po, such that 
Pi = Cipo and p r = C r pp, a dim ensionless form of the 
current is given by y/2-!rm/ (k B T)j/ po. This is plotted in 
Figure |3 for the cases where the reservoir densities are 
equal (C/ = C r ) or different (C; < C r ). Figure [5] shows 
that the current saturates for large \a\. For positive a, 
the current at saturation is due exclusively to ions from 
the left reservoir and y / 2irm/(k B T)j / p — ► Ci\ for nega- 
tive a, ^2irm/(k B T)j / po — > C r . For values of \a\ below 
saturation, the magnitude of the current increases as the 
density of traps increases. This is because traps reduce 
the negative current contribution j r from the right reser- 
voir (for a > 0), without affecting the current ji of ions 
moving from the left, as can be seen in the insets, where 
j r and ji (in dimensionless form) are plotted individu- 
ally for the case where pL = 1. The discontinuity in 
the current at a = 0, observed for finite concentrations 
of traps (p > 0), reflects the fact that the model is no 
longer valid in the absence of a field, when there is no 
stationary solution (since ions that are stopped by a trap 
are not then re-accelerated). In the case of unequal ion 
densities in the two reservoirs, the current-voltage curves 
are asymmetric, as shown in Figure ^jp. The saturation 
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FIG. 2: Dimensionless current -J '2irm / '(fesT)j / po as a func- 
tion of maL /(kBT), where the reservoir densities are pi = 
Cipo and p r = CVpo, for various values of the dimensionless 
stopping trap density pL. Solid lines: pL — 0, dotted lines: 
pL — 1, dashed lines: pL — 2. (a): Equal reservoir densi- 
ties, p r = pi; C r = 1, Ci = 1 (b): p r = 2pi; C r = 2, Ci = 1. 
The insets show the currents ji and j r (in dimensionless form) 
due to ions originating in the left (circles) and right (squares) 
reservoirs, for the case pL — 1. Values for the current in ab- 
solute units can be obtained by substituting absolute values 
for the physical parameters a, m, L, p r , pi, p and ksT . 



value of | j | is now larger for negative a, and j is negative 
for small positive values of a. 

Thus far, we have assumed that any ion which encoun- 
ters a trap is stopped, regardless of its velocity. How- 
ever, ions with low kinetic energy could be expected to 
be more likely to be bound by a polar residue lining a 
nanopore, than those with more energy. We now con- 
sider a variation on our previous model, in which a sin- 
gle trap is present at position x = xq inside the pore 
(—L/2 < xo < L/2), which presents an "energy barrier" 
of height Eq — mv^/2 to all ions crossing x — xq. Wc 
shall consider two possible modes of action of this trap. 

In model A, the trap at x — xq stops all ions with 
kinetic energy below the barrier height: mv 2 /2 < Eq 
(and subsequently releases them to be re-accelerated by 
the electric field), but has no effect on ions with energy 
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mv 2 /2 > Eq. The appropriate kinetic equation reads: 



S(v) / dw\w\f(x,w) - 0(v o - \v\)\v\f(x, v) 



(16) 



Equation (flfcip) can be solved analytically, but we shall 
instead use simple arguments, as before, to obtain the 
current j without the explicit form of f(x, v). As above, 
any ion entering the channel at the left extremity (x = 
—L/2) with velocity v > will eventually reach the right 
extremity (x = L/2), regardless of whether it is stopped 
by the trap. The contribution ji of these ions to the 
current is therefore given by l|10|) . However, ions entering 
the channel at x = L/2 with velocity v < will only reach 
x = —L/2 (and hence contribute to the current) if they 
are not stopped either by the field or by the energetic 
trap at xq. There are two possibilities, depending on the 
barrier height Eq: 

(1) Eq > ma(xQ + L/2). In this case, any ion which 
reaches the trap with energy greater than Eq, and so 
is not stopped, must have sufficient energy to overcome 
the remaining part of the opposing field between Xq and 
—L/2. In order to have energy > Eq on reaching the 
trap at xq, an ion must enter the channel at x = L/2 
with velocity v such that: 



mv 



>E 



Xq 



(17) 



so that, assuming an incoming Maxwell distribution, the 
distribution function of these ions is: 

f r (x, v) = p r 9(-v) 6 (y - ax - - ax a ) \ (18) 



2nk B T 



exp 



rn 



v 
~2 



Using © we find that in this case the contribution to the 
current due to ions from the right hand reservoir is: 



Jr 



2nm 



exp 



Eq 



k B T k B T \ 2 



Xq 



(19) 



(2) Eq < ma (xq + L/2). In this case, ions which pass 
through the barrier at xq do not necessarily have suffi- 
cient energy to overcome the remaining part of the field 
between xq and —L/2. The only ions which contribute to 
the current are those entering the channel with velocity 
V, such that: 



mv 



> maL 



(20) 



These make a contribution to the distribution function: 



fr(x,v) = p r 9(-v)t 



aL (21) 



2itk B T 



exp 



m 
~k~RT 



V v 2 



which results in a contribution to the current: 

maL 1 



Jr 



-Pr 



k B T 



■ exp 



(22) 
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The total current j = ji + j r for model A is therefore: 
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Eq — ma 
Eq 



L 
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(23) 
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L 



-p r 9 ( ma ( xq + — ) - Eq 



ma 
k^f 

exp 



2 ~ Xo 



maL "I 

JbY j 

Expression (|23|l is, of course, only valid for positive values 
of a. The equivalent expression when a < can easily be 
shown to be: 



Ja<0 



k B T 
2irm 



p r + piO E 



x exp 



■piH \ -ma [ - 



Eq 



k R T k n T 



Xq 



ma\ — - x 



L 



(24) 



2 +XQ 



E exp 



maL 1 



The current-voltage curves for model A are shown in 
Figure|31 in dimensionless form as in Figure|2 For clarity, 
we consider the case where only the right-hand reservoir 
contains ions: p\ = 0. In Figure the position of the 
trap is fixed at xq — and the barrier height Eo/k B T 
is increased. When maL/k B T exceeds the critical value, 
given by maL = E /(l/2 + x /L), j no longer depends 
on Eq and all the curves become identical. However, in 
the regime where maL < Eq/(1/2 + Xo/L), j depends 
strongly on Eq, being increased on increasing the bar- 
rier height. This can be easily understood, since ions 
which come from the right-hand reservoir and are im- 
peded by the barrier make a negative contribution to 
the current. An interesting general observation can be 
made here, that the presence of an energetic barrier can 
have the effect of increasing the ionic current. In Fig- 
ure Ob, the barrier height is fixed (E /k B T = 1) and 
the trap is moved towards the right-hand end of the 
pore. The current-voltage characteristics are seen to be 
rather sensitive to the position of the trap in the regime 
maL < Eq/{\/2 + xq/L), although there is no depen- 
dence for larger maL/k B T. As xq increases, j decreases, 
although the value as a — > remains unchanged. 

We also consider an alternative energetic barrier 
model, model B. Here, ions encountering the trap with 
energy greater than the barrier height, mv 2 /2 > Eq, do 
not continue unperturbed, as in model A, but instead 
lose energy Eq, being relea sed by t he trap with reduced 
velocity v', where \v'\ — \Jv 2 — Vq. Less energetic ions 
with mv 2 /2 < Eq are stopped by the trap, as in model 
A. Following a line of reasoning as for model A, one finds 
that the total current (when a > 0) is: 



J 



k B T 
2nm 



pi - p r exp 



k B T 



aL + ^ 
m , 



(25) 
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2 
maL/ksT 



FIG. 3: Dimensionless current \J 2Tvm/(kBT)j/po for Model 
A (Equations 1231 and 1241 ') as a function of maL/(kBT), (po 
defined as in Figure for p r = po; pi = 0. (a): Energetic 
trap fixed at centre of pore, xo/L = 0. Solid line: _Bo/(fcflT) = 

0, dotted line: E /(k B T) = 0.5, dashed line: E /(k B T) = 

1. (b): Height of barrier fixed, E /{k B T) = 1. Solid line: 
xo/L = 0, dotted line: xq/L — 0.2, dashed line: xq/L = 0.4. 



Note that for model B, the current j does not depend on 
the position xq of the trap. Comparing expression H25|) 
with H14fl and l|15|) . we see that the current in model B 
with energy barrier Eq is identical to that through the 
channel with stopping traps investigated at the beginning 

of this section, if E /k B T = j_^, 2 p( x ') dx' . 



IV. FINITE CHANNEL WITH FRICTION 



Equation (|26|l may be solved subject to boundary condi- 
tions specifying the incoming particle fluxes from the left 
(x = —L/2) and from the right (x = L/2), i.e.: 



rf [ -,v ) ilr = - ■/>, 



k B T 
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(27a) 



(27b) 



In the limit of an infinitely long channel (L — » oo), 
the distribution must be homogeneous (f(x,v) — > f(v)); 
the solution of the corresponding Fokker-Planck equation 
(i.e. <j2"fj)l without the vd/dx operator in the free flow 
term) is: 



exp ■ 



mi a 

2k B r y~ 7 



(28) 



On the other hand, a particular inhomogeneous solution 
of Equation l|26|) in a finite channel is: 



f{x,v) 



exp 



(29) 



We now look for a solution of the Fokker-Planck equation 
(|26|l for finite L, satisfying the boundary conditions 12711 . 
in the form of a linear combination of the two solutions 
BB> and (H3): 
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k B T 
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(30) 



Distribution l|3U|) indeed satisfies l|26|l for all values of the 
coefficients A and B. Imposing the boundary conditions 
()27J) . we obtain: 



£ = 



X 



(31) 



We now turn to a model where no stopping traps are 
present (p{x) — 0), but ions undergo frictional collisions 
inside the pore (of finite length L). The thermalizing 
effect on the ion of these collisions with the channel sur- 
face and with other molecules (e.g. water) is modelled by 
a Langevin mechanism, represented by a Fokker-Planck 
operator, so that the stationary state kinetic equation 
now becomes: 



' d d \ d 

v ox ov I ov 



k B T d 
m dv 



f(x,v) 
(26) 



where 



_ . . / maL \ 
X = 2 smh I — — — ] exp 



\2k B T / 



a 2irm 
+ 7 V fc^T 

and 

A = 



cosh 



1 



maL 
2k B T 



2 sinh ( 



2k B T 



2k B T~f 2 

( maL 



sinh 



\2k B T 



Pr - Pl + 



a / 27r?n 



7 V k B T 



B 



(32) 



(j> (v) dv 



(33) 
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The current through the channel can then be calculated 
using Equation (JSJ: 



J = -B 

7 



(34) 



Expression (|34|1 for the current simplifies greatly in the 
limit of vanishing applied field (a — > 0), when the only 
driving force is diffusion under the action of the density 
gradient (pi — Pr)/L. Substituting 1)31(1 into (|34Jl ■ one 
finds: 



1- ■- / k EiT pi - p r 

a™ 3 ~ V 27rm (1 + 7L0 T (O)) 



(35) 



i.e. the friction reduces the current by a factor 
1/ (l + r yL(j) T (0)) compared to the free flow result J7J). 
The a — > limit of the distribution function f(x,v) a ^o 
and the resulting density profile are discussed in Ap- 
pendix [5] f(x,v) a ^o for a given L and 7 turns out not 
to be everywhere positive, pointing to a fundamental dif- 
ficulty in applying the Fokker-Planck equation (|26[) in a 
system of finite spatial extension L. 

Another interesting case is the limit of strong friction 
(7 — » 00). B then becomes: 



lim B 

7— >oo 



pi - p r exp 



J maL \ 

\ k B T j 



(36) 



and the resulting current, given by inserting l|36l) in (1341) . 
reduces to the classic GHK expression |23j, |24|, which 
arises from solving the 1-d diffusion equation in a con- 
stant external field |6|]. 

GHK theory predicts that the current increases lin- 
early with voltage across the channel for large voltages. 
However, the behaviour of j for large a in this model 
is considerably different: fixing 7 and taking the limit 
a — > 00 in Equation l|31Jl . we find that the current satu- 
rates for large applied fields: 



lim j 

a — >oo 

lim j 

a — > — 00 



Pi 



k B T 
27rm 



k B T 
2irm 



(37a) 
(37b) 



Equations H37fl . which are identical to the saturation val- 
ues of the current for the models presented in section 
lllll correspond to the situation where all ions crossing 
the channel in the direction of the field contribute to the 
current and all ions attempting to penetrate the channel 
against the field are turned back. 

An important experimental quantity is the "reversal 
potential" : the voltage across the channel for which the 
total ionic current is zero. In the case where the ionic 
species are the same in the two reservoirs, the accelera- 
tion ao at which the current is zero is given by cancelling 
the numerator of (|3ip. which yields: 



ao 



k B T 
mL 



hi 



(38) 




-2 2 

maL/ksT 



FIG 

tion 



4: Dimensionless current ^/27rm/(fc_BT)j/po as a func- 
of maL/ '(ksT), for values of dimensionless friction 
L-y^/m/(k B T) of 0.1 (solid lines), 1.0 (dotted lines) and 10.0 
(dashed lines), (a): Equal reservoir densities p r = pi = po; 
inset shows results predicted by GHK theory, (b) Reservoir 
densities p r — po; pi = (above) and p r = 0; pi = po (below); 
again, inset shows results of GHK theory. 



Expression l|38() is identical to the GHK prediction. How- 
ever, if the reservoirs contain different species: for exam- 
ple potassium on one side and sodium on the other, the 
model will no longer agree with GHK theory. 

Plots of the dimensionless current ^J2irm/ (k B T)j / p^ 
versus maL / (k B T) are shown in Figure for values of 
the dimensionless combination Ljy^m/ (k B T) of 0.1, 1.0 
and 10.0. In both Figure 0^, and Figure 0Jd, the insets 
show the results of GHK theory. Figurc^Ji shows the cur- 
rent through the channel when the ionic concentrations 
in the two reservoir are equal (p r — pi). While GHK 
theory predicts linear asymptotic behaviour, the current 
given by Equation l|34|) shows saturation as \a\ — > 00. As 
the friction coefficient decreases, the current-voltage rela- 
tion becomes steeper and deviates further from the GHK 
results. In Figure 0Jd, the current shown in Figure^ is 
divided into the contributions of ions originating in the 
right (shown above) and left (shown below) reservoirs. 
As expected, for large positive a, the current is due only 
to ions from the left, and for large negative a, it consists 
only of ions from the right. 
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V. INFINITE CHANNEL WITH FRICTION 
AND TRAPS 

We next address the full version of the model system 
described in section [H] single ions moving under the in- 
fluence of constant accelerating field, a Fokker-Planck 
thermalizing mechanism and stopping traps. We shall 
consider only the case where the velocity of an ion en- 
countering a trap is set to zero, irrespective of its initial 
velocity, and where the average distribution of the traps 
is uniform (p(x) — p). An analytic solution of the ki- 
netic equation is presented in the limit of an infinitely 
long pore [L — > oo). This solution may prove useful in 
analysing ion flow through carbon nanotubes or the long 
pores found in zeolites. 

We first introduce dimensionless position and velocity 
variables y and u: 



1 

-y 

P 



k B T 



(39a) 
(39b) 



as well as a dimensionless distribution function F(y, u) 
through the transformation 



f(x,v)dxdv = f l-y, 
\ P 




dy du 



= F(y,u)dydu (40) 
Using H39fl and (|40p . the kinetic equation (JTJ becomes: 



P { S(u) 

d_ 

du 



+oc 



dw\w\F(y, w) 
F(y,u) 



\u\F(y,u) 



0_ 

du 



where the dimensionless coefficients a and (3 are defined 
by: 



a 







kef 



(42a) 



(42b) 



We were unable to solve the inhomogeneous equation 
(|41|l analytically. An analytic solution may, however, be 
obtained in the in the limit of an infinitely long pore 
(L — > oo), when the ion distribution no longer depends 
on y and the problem is spatially homogeneous. Equation 
(|41|l then simplifies to: 



dF(u) 
du 



P\S(u) 



+ 00 



dw\w\F(w) — uF(u) 



du \ du ' 



(43) 




FIG. 5: Current J = m/ (kBT)j as a function of a, for 
values of f3 of 0.0 (solid line), 0.25 (dotted line), 0.5 (dashed 
line) and 0.75 (dot-dashed line). Inset: Total J (dashed line) 
as well as components of J towards the right (circles) and 
towards the left (squares), for the case where /3 = 0.5. 



The solution of Equation l|43|l is obtained as sketched in 
Appendix IO The result is: 



F(u) = Aexp 



(u — a)' 



(44) 



9{u) D 0{l3+a) (a+2f3) D 0(0 _ a) (u-a+2f3) 
+6(-u) (-a+2/3) D w+a) (~u+a+2f3) 



(41) where the D p (z) are parabolic cylinder functions. The 
constant A determines the number density of ions inside 
the channel (in a finite channel this would be set by the 
reservoir densities): here, we assume one ion per unit 
channel length, so that A can be obtained numerically 
from the normalization condition f^°° F(u)du = 1. 

In the limit (3 — > 0, i.e. in the absence of traps, 
|@2J reduces to the result ((2HJ (noting that D (z) = 
exp [— 2 /4] and reverting to dimensional units), and the 
current is linear in a. In the limit a — > 0, i.e. in the 
absence of acceleration, the solution (|44|l is seen to be 
an even function of u, so that the current j vanishes, 
as expected for an infinitely long, spatially homogeneous 
channel. 

In the general case, when both a and (3 ^ 0, the 
current must be calculated by numerical integration, af- 
ter substituting (|44() in (|5j|. Figure |5] shows the current 
J = uF(u)du = \Jm/ (fcsT)j', as a function of a 
for (3 between and 0.75. The inset shows the forward 
and backward components of J when f3 = 0.5 (given 
by integrating over the coefficients of 9(u) and 0(—u) in 
Equation (03J). There is a qualitative difference in the 
behaviour of the current when f3 = 0, where the relation 
between J and a is linear, as noted above, and when 
(3 7^ 0, where it is non-linear. Thus even a very small 
density of stopping traps (for example, due to defects 
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or impurities) can have a dramatic effect on the current 
flowing through the pore. On estimating typical values 
of the physical parameters a, 7, m and p, we find that 
a and (3 are in fact likely to be small, perhaps of order 
0.01-0.1. 



VI. FINITE CHANNEL WITH THERM ALIZING 
CENTRES AND TRAPS 

In Sections IIVI and ^ the effect of friction and ther- 
malization on the motion of single ions was modelled by 
the Fokker-Planck collision operator. For the channel 
of finite length, this leads to the fundamental problem 
that imposing the incoming ion fluxes from the reser- 
voirs at both ends of the channel results in a stationary 
distribution function f(x, v) which is not positive definite 
(see Appendix [BJ). In this section, we therefore replace 
the Fokker-Planck mechanism by an alternative thermal- 
ization process. We consider a model in which the 1-d 
channel contains N "thermalizing centres" , at positions 
Xi, such that: 



-— < X\ < x 2 < 



< x N < — 



(45) 



When an ion reaches a thermalizing centre, its incoming 
velocity v is replaced by a new velocity v' drawn from a 
Maxwell distribution <p T (v'), Equation ©. The channel 
also contains a series of "energy barriers", of the type 
denoted "Model A" in section ITTT1 a barrier of height Ei 
temporarily stops an ion with kinetic energy mv 2 /2 < Ei 
but has no effect if mv 2 /2 > Ei. An energy barrier of 
height Ei is located between each pair of neighbouring 
thermalizing centres at x^i and X4: if Ei = 0, this is 
equivalent to having no energy barrier present. In be- 
tween encounters with thermalizing centres and energy 
barriers, ions move with constant acceleration a, which 
is taken to act towards the right. 

We now analyze the stochastic process defined by this 
model, leading to an exact calculation of the stationary 
current j. The key quantity is the probability p(i — > i + 1) 
that an ion which is thermalized by the centre at Xi, 
subsequently encounters the next thermalizing centre at 
Xi + \ {i.e. it reaches Xi+i before Xi-x). We first note that 
an ion which leaves the centre at Xi, moving towards the 
left, requires a minimal energy e(i,i — l) to penetrate the 
field and energetic barrier and reach the centre at Xi-±, 
where: 



e(i, i — 1) = ma{xi 



(46) 



We shall adopt the convenient notation xo = —L/2 and 
xn + i = L/2, so that Equation i|4tj|) remains valid for 
i = l and for i = N + 1. Ions leaving Xi towards the left 
with energy less than e(i, i — 1) will be stopped and re- 
accelerated towards the right, returning to Xi. Since ions 
leave the thermalizing centre with a Maxwell velocity dis- 
tribution, the probability w(i,i — 1), that an ion leaving 



Xi (in either direction) has energy less than e(i, i — 1) is 
given by: 



2 j.^e(Li-l)/k B T 

w(i,i-l) = —= duexp(-u 2 ) (47) 

V 71 " Jo 

An ion that leaves Xi towards the right, on the other 
hand, will certainly reach the thermalizing centre at ij+i, 
regardless of whether it is stopped by the energetic bar- 
rier Ei + \. Thus on leaving Xi, an ion may be sent to the 
right (with probability 1/2), and reach x$+i, or it may 
be sent to the left, be stopped and return to Xi (with 
probability w(i, i — l)/2), or lastly it may be sent to the 
left and reach Xi_\. 

The probability p(i — ► i + 1) that an ion leaving n 
reaches Xi+i before Xi-x can be found by summing over 
all the possible ways that this might happen. The n-th 
term in the series corresponds to the scenario where an 
ion is sent n times to the left (i.e. towards Xi-i) and 
returns to Xi before eventually being sent to the right 
(i.e. towards Xi+x). We thus obtain a geometric series: 



p(i 



i + i) = i + i 
' 2 2 



w(i, i — 1) 


1 


w(i, i — 1) 


2 


+ 2 


2 



2 - w(i,i- 1) 



(48) 



The probability for the transition in the opposite direc- 
tion (xi to Xi-x) is then clearly: 

p(i-+i-i) = i-p(i^i + i) = 1 ~ w f: i r]\ (49) 

2 — w(i, i—l) 

The probabilities l|48|l and 14911 may now be used to de- 
termine the probability % that an ion starting from x^ 
eventually leaves the channel through the right end at 
xjv + i = L/2. The qi satisfy the (detailed balance) equa- 
tions: 



q% =p(i -> i- + p(i -* 

with the boundary conditions 

<7o = 0, q N +i = 1 



1)%- 



(50) 



(51) 



Defining the differences = — ft-i, one finds from 
(|5Ujl that: 



A 



i+l 



p(i 



1) 



1 = 



p(i 



1) 



p(i 



1) 



(52) 



Taking the product of both sides of Equation (|52|l over 
1 < i < n leads to: 



n 

i=l 



A H 



1n+l 



n 



p(i 



1) 



^ p(i -> i + 1) 



(53) 



Summing both sides of the second equality in Equation 
(|53|l over 1 < n < N, we arrive at: 



1-gi 

9i 



N n 



n=l i=l ^ y 



i-D 



i + l) 



(54) 
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The only ions which make a contribution to the current 
are those which come from the left reservoir, pass through 
the whole channel and exit at the right end, and those 
which come from the right reservoir and exit at the left 
end. We now calculate the probability p(—L/2 — ► L/2) 
that an ion entering the channel from the left reservoir 
will exit through the right end, and thus contribute to 
the current. On entering the pore at —L/2, the ion will 
reach the thermalizing centre at x± with probability 1, so 
that the definition of qi, together with Equations 1)540. 
gHJ and dJJ}, lead to: 



JV n 



p(-L/2-»L/2) = = i + snsj 



p(i — * i — 1) 



^— ; 11 p(i — * i + 1) 



(55) 



{JV n V 
i + EIIt 1 -^^" 1 )] 
n=li=l J 



Consider next an ion entering the channel from the right. 
It will reach the thermalizing centre at xn with proba- 
bility exp[-e(7V + 1, N)/k B T}. Therefore: 

p(L/2 -» -L/2) = {l-q N )exp[-e(N+l, N)/k B T] (56) 

We find qN by setting n = N in Equation (|53|) and using 
Equations Jjjj and (jSTIj) : 



1 - ?JV 



9i 



n 



p(i — ► i — 1) 



11 p(i 

= p(-L/2-> L/2) JJ[1 -«;(*, 



N 



(57) 



i=l 



We now combine Equations (|56[1 . I|56|) and l|57|l and con- 
clude that ions coming from the right reservoir contribute 
to the current with probability: 



p(L/2 -> -L/2) = 



nL[l--(M-l)]exp{^^} 



w(i,i- 1)] 

(58) 

The stationary current is given by the sum of the incom- 
ing fluxes from the left and right reservoirs (found by as- 
suming incoming Maxwell distributions), multiplied by 
the probabilities p(-L/2 -> L/2) and p(L/2 -> -L/2), 
which determine the extent to which the incoming flux is 
reduced by the action of the thermalizing centers: 

3 = 31 + 3r (59) 

[ Pl p(-L/2 -> L/2) - Pr p(L/2 -> -L/2)] 



fcjjL 

27TTO 



Inserting Equations (|56() and l|58|) into ijSSJ, we find: 



k B T 
2ixm 



(60) 



Pj - Pr IltJ 1 - * - l)]exp[-e(iV + 1, N)/k B T] 
1 + Eli niLi [1 -«;(»,«-!)] 



Combined with formulae (|46(l and (|47|l . Equation i|61|) 
provides an explicit expression for j as a function of pa- 
rameters defining the internal structure of the channel. 

In the absence of applied field (a — 0) and energy 
barriers (all Ei = 0), w(i, i — 1) = and the current l|61|l . 
now due to the effect of the thermalizing centres only, 
takes the particularly simple form: 



k B T 
2nm 



Pi - Pr 



N+l 



(61) 



i.e. both incoming fluxes are reduced by the same factor 
l/(iV + 1). The prediction i|35|) of the FP equation under 
the same conditions (a — Ei — 0) coincides with 161II . 
provided: 



— = ■y\/m/2Trk B T 

1j 



(62) 



The "effective friction" introduced by the thermalizing 
centres is thus proportional to their density. Physically, 
()62J) also means that the relaxation time 7 -1 is of the 
order of the time taken by an ion to cover the average 
distance L/N between the thermalizing centres with ve- 
locity yJk B T /m. This equivalence between the FP and 
thermalizing centre results does not hold, however, in the 
presence of an accelerating field (a / 0). 

As in the case described by the FP collision term (cf. 
Equation lpF7)ll. the current j saturates for large applied 
fields (a — > 00) at the value ji. This is because all ions 
coming from the left are driven through the channel by 
the strong field, while no ions are able to cross the channel 
successfully from the right. 

Formula (|61() simplifies greatly when the thermalizing 
centres are evenly distributed and in the absence of en- 
ergy barriers, i.e. when (xi — 2^-1) = L/(N + 1) and 
Ei = for all 1 < i < N. In that case: 



1 - s 



N 



k B T 
2-Kvn 



Pi 



N 

p r s exp 



maL 



(N+l)k B T 



where 



_ du exp(— u 2 ) 



(63) 



(64) 



Figure shows the effect on the dimensionless current 
^J2i:m/ {k B T)j / pq of increasing the number N of evenly 
spaced thermalizing centres, when there are no energetic 
barriers {Ei = for all i). In Figure^,, the reservoir 
densities are equal, pi — p r = Po- As N increases, the 
current decreases, requiring larger values of maL / (k B T) 
to approach its asymptotic value. The inset shows the 
contributions to the (dimensionless) current from the left 
and right reservoirs when N — 5. Figure shows re- 
sults for the same channel, when the density of ions in 
the right-hand reservoir is twice as large as that in the 
left: p r = 2pi = 2p . In this case, the negative contri- 
bution j r is increased and dominates for small values of 
maL/{k B T). 
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FIG. 6: Dimensionless current ^/ 2nm / (ksT)j / po as a func- 
tion of maL/(kBT), where the reservoir densities are pi = 
Cipo and p r = C r po, for channels containing an increasing 
number N of evenly spaced thermalizing centres. Energetic 
barrier heights Ei are all set to zero. Solid lines: N = 0, 
dotted lines: N — 1, dashed lines: N = 5, dot-dashed lines: 
N = 10. (a): Equal reservoir densities, p r = pi; CV = 1, 
Ci = 1 (b): p r = 2p;; CV = 2, C; = 1. The insets show the 
currents ji and j r (in dimensionless form) due to ions origi- 
nating in the left (circles) and right (squares) reservoirs, when 
N = 5. 




°0 2 4 6 8 10 

(b) maL/k B T 



FIG. 7: Dimensionless current v/ 27rm / (&b T) j / po as a func- 
tion of maL/(fesT), for equal reservoir densities pi = p r — po- 
(a): Channel contains N — 5 thermalizing centres, evenly 
spaced in the range —6/2 < x < b/2. All barrier heights Ei 
are set to zero. Solid line: b — L, dotted line: b — L/2, dashed 
line: b = L/4. (b): Channel contains N = 4 thermalizing cen- 
tres, evenly spaced in the range —L/2 < x < L/2, and one 
energetic barrier E3 > 0, located between the central pair of 
thermalizing centres. Solid line: Ez/{ksT) = 0, dotted line: 
E 3 /(k B T) = 1, dashed line: E 3 /(k B T) = 2. 



We have also investigated the effect of changing the 
spatial arrangement of the thermalizing centres, once 
again in the absence of energetic barriers (Ei = for 
all i). In Figure [7Ji, the channel contains N = 5 
thermalizing centres which are all located in the range 
—b/2 < x < b/2, where b < L. Within this range the 
thermalizers are evenly spaced. Results are shown for 
equal reservoir densities, pi = p r = po- As b decreases 
and the thermalizers become more localized in the middle 
of the pore, the current increases, approaching its asymp- 
totic value for smaller values of maL / '(ksT). However, 
the results for b = 0.5L (dotted lines) and b = 0.25L 
(dashed lines) are rather similar, indicating a limiting 
current-voltage relationship for small b. 

Interesting effects are obtained on including energetic 
barriers. Figure[7|3 shows results for a channel containing 
= 4 thermalizing centres, with a single barrier Ei, lo- 
cated between the second and third thermalizers (i = 3), 
i.e. in the central one of the 5 possible positions. Once 



again, pi = p r = po- As the barrier height E3 is in- 
creased, the current increases, showing that inserting an 
impedance to ion passage can actually enhance the total 
ion flow through the channel. This apparently somewhat 
counter-intuitive result can in fact be easily understood. 
Let us consider an ion that is released by the thermal- 
izing centre at x\. If it is sent out towards the right, 
the ion will certainly reach a^+i, while if it is sent to the 
left, it will be stopped and return to Xi with probability 
w(i,i — 1), which is an increasing function of Ei. Thus 
increasing Ei increases the chances of an ion eventually 
arriving at Xj+i, and thus enhances the current. This 
phenonenon may be of interest for biological ion chan- 
nels, where the selectivity filter might play the role of an 
energetic barrier. 
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VII. CONCLUSION 

In this paper, we have introduced some simple kinetic 
models for the transport of independent ions through 
narrow pores, under the influence of a constant acceler- 
ation, due to an applied external electric field. We con- 
sider only one-dimensional motion along the pore axis. 
Analytic results have been obtained for the stationary 
ion current j and, in some cases, also for the stationary 
non-equilibrium distribution function f(x,v). Boundary 
conditions at the ends of the pore are consistent with 
the presence of reservoirs containing ions at equilibrium, 
which determine the flux of ions entering the pore. The 
models include traps or energy barriers, which represent 
the temporary binding of an ion to polar residues lin- 
ing the pore surface, and they also account for the fric- 
tion and thermalization due to collisions of the ions with 
molecules {e.g. water) inside the pore as well as with its 
inner, confining surface. Initially, this friction and ther- 
malization was included using a Fokker-Planck operator 
in the kinetic equation. We were unable to find a solu- 
tion for a pore of finite length L when both traps and the 
Fokker-Planck description of friction are present. How- 
ever, we have derived an explicit solution for the homo- 
geneous case of an infinitely long pore. For pores of finite 
length, solutions are given for models with stopping traps 
or energy barriers in the absence of friction. For stopping 
traps, the current does not depend on the spatial distri- 
bution of the traps, but for a single energy barrier, there 
may be a dependence on its position, depending on its 
mechanism of action. We have also solved the kinetic 
equation for the distribution function f(x,v), for a finite 
channel in the case where the Fokker-Planck mechanism 
is present but there are no traps or energy barriers. In 
this case, on imposing the incoming fluxes from the reser- 
voirs at both ends of the pore, f(x, v) turns out not to 
be positive definite for short channels and/or small val- 
ues of the friction constant 7. This unphysical behaviour 
can be understood in terms of competing time scales for 
high velocity ions (which pass through the channel before 
they can be thermalized) : the resulting ion current re- 
mains well-behaved, as does the number density profile. 
In view of this deficiency of the Fokker-Planck mecha- 
nism in a pore of finite length, we have introduced an 
alternative model, whereby ions are instead thermalized 
by encounters with a series of "thermalizing centres" , lo- 
cated at given positions inside the pore. Energy barriers 
may also be present. 

An important conclusion arising from all the models 
which were considered is that the current j invariably 
saturates as a function of the external field (or equiv- 
alently the constant acceleration a), since it is limited 
by the incoming flux from the reservoirs. This satura- 
tion behaviour contrasts with the linear increase of j with 
voltage predicted by the classic GHK result, which can be 
derived by solving the one-dimensional diffusion equation 
in the presence of a constant external field. A further in- 
teresting observation that emerges from this work is that 



the presence of stopping traps or energy barriers inside 
the pore can increase rather than decrease the ionic cur- 
rent. This is because the steady state flow of ions crossing 
the pore in the direction of the applied acceleration is un- 
affected by the traps or barriers, while the flow of ions 
against the field is reduced. When the former contribu- 
tion is the dominant one, the current will be enhanced 
by the traps or barriers. 

A question that arises is whether there is any corre- 
spondance between the two models of friction and ther- 
malization considered in this work: the Fokker-Planck 
mechansim and the the model involving N thermalizing 
centres. In the absence of an applied field and of energy 
barriers, when ionic motion is driven only by the concen- 
tration gradient across the pore, an equivalence can be 
established between the two models (c.f. Equation l|62|)l. 
However, under more general conditions, we have found 
no one-to-one correspondance between the two descrip- 
tions of dissipation. 

Future work will include a complete numerical analysis 
of the model involving both energy barriers and thermal- 
izing centres. We also plan to extend the kinetic mod- 
els to include the possibility of collective ion permeation 
through a pore, by including ion-ion interactions. These 
are believed to play an important role in ion transport 
through some biological channels 0, 0, HTl IT^ . Appro- 
priate selection of parameters such as the pore length L, 
the friction coefficient 7 or the number and positions of 
the thermalizing centres and the positions and heights 
of the energy barriers, to correspond to the structure of 
real ion channels, should allow the predictions of these ki- 
netic models to be compared to measured current-voltage 
characteristics. 

APPENDIX A 

Here we derive an analytic solution for the kinetic 
equation @, for the case where p{x) = p: i.e. we find the 
distribution function f(x,v) for ions experiencing a uni- 
form external electric field and traps distributed accord- 
ing to a Poisson law corresponding to a uniform average 
density p. 

The distribution function f(x,v) is split into the con- 
tributions f + (x,v) and f~(x,v) of ions moving to the 
right and to the left: 

f(x,v) = 6(v)f + (x,v) +6(-v)f-(x,v) (Al) 

Substituting 1A1() into © (with p(x) — p), we obtain: 

^v^- + a^- + pv^j f+(x,v)=0 (A2a) 

{ v iL +a ^~ pv ) r{x,v) =o (A2b) 

Equations (IA2|) imply that: 

f + {x,v) = exp[-px]F+ (v 2 /2-ax) (A3a) 
f-(x,v) = exp[px]F~ (v 2 /2~ax) (A3b) 
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The as yet unknown functions F + and F~ are linked by 
the requirement JSJ) that the current j be independent of 
x: 



j = I dvvf(x,v) — constant 



(A4) 



Substituting l|A3jl into 1|A4J) . and defining w = v 2 /2 — ax, 
we obtain: 



j = constant 



(A5) 



/oc />oo 
F + (w) dw — exp [px] / F~(w) 
-ax J —ax 



dw 



Multiplying l|A5J) by exp [px] and differentiating with re- 
spect to x, we obtain: 



aF + (—ax) — jpexp [px] 

t>oo 

+ 2p I F~(w)dw + aF~(-ax) 



(A6) 
exp [2px] 



Equation i|A6|l is valid inside the channel, i.e. for values 
of x in the range —L/2 < x < L/2. The argument —ax of 
F + and F~ therefore ranges between —aL/2 and aL/2, 
so that the relationship (|A6I) between F + (w) and F~(w) 
holds for —La/2 < w < La/2. Since w = v 2 /2 — ax, this 
corresponds to v 2 /2 < a{L/2 + x). 

The distribution of particles moving against the field, 
f~(x,v), can be obtained using simple arguments. If a 
particle reaches position x with negative velocity v, then 
it must have had energy mu 2 /2 — mv 2 /2 + ma(L/2 — x) 
at the right-hand channel entrance, and have encountered 
no traps over the distance (L/2 — x), since on encoun- 
tering a trap an ion is re-accelerated by the field towards 
the right. Assuming a Maxwell distribution for the ve- 
locities of the incoming ions and noting that the Poisson 
probability for encoutering no traps over this distance is 
exp [— p(L/2 — x)], we find that: 



/ (x,v) = e^cp[px]p r 



2nk B T 



(A7) 



exp 



k B T 

and thus from (IA3I 

F-(w)= Pr 



(v 2 /2 . 







ma \ 


exp 




k B Tj\ 



2irk B T 



exp 



exp 



(A8) 

Substituting the result l|A8|) into (| A6|> . we find that for 
v 2 /2 < a(L/2 + x): 



■ P 

= J— exp 
a 


P 2 

2a 


Pr_ 




m 


2p- 




T 


a 


k B T 


a _ 


V 2itm 


x exp 




m ( if 


\ 

— ax 


\ 


exp 


px - 


P 2 
V 

a 






k B T \ 2 







x exp 



2{ P+ k^ 



(A9) 



Using a similar argument to that above, we can determine 
the remaining part of f + (x,v), for v 2 /2 > a(L/2 + x). 
Ions with v 2 /2 > a(L/2 + x), moving towards the right, 
cannot have been stopped by a trap, since their en- 
ergy is larger than that due to the accelerating field 
over the distance travelled in the pore. These ions 
must have entered the pore at x = —L/2 with energy 
mu 2 /2 = mv 2 /2 — ma(L/2 + x), and have encountered 
no traps over a distance L/2 + x. Assuming a Maxwell 
distribution of incoming particles, we obtain a contribu- 
tion to f + (x, v) of: 

9 (v 2 /2 - a (L/2 + x)) exp [-p (L/2 + x)] x (A10) 
(v 2 /2- a(L/2 + x)) 



Pi 



in 



2nk B T 



exp 



kuT 



Noting that at x = —L/2, v 2 /2 > a(L/2 + x), we can use 
equations l|A7|) and ijAlOfl for x = —L/2 to calculate the 
current: 



J 



k B T 
2itm 



pi - p r exp 



—L 



kuT 



(All) 



Substituting this result in (|A9 
distribution function f(x,v) is: 

f(x,v) = 6(-v)p r (j) T (v)exp 



the final result for the 



+ P (a -L/2) 



k B T 



+6(v){9 (v 2 /2 - a (L/2 + x)) Y 
+9 (a (L/2 + x)-v 2 /2)Z) 



(A12) 



where 



Y = pi(f> T (v) exp 
Z = 



ma 
k^f 



k R T 



2nm 



pi - p r exp 



-pUx + L/2) 

( ma 
L ( P+ k^T 



(Al3a) 



P 

x - exp 

a 

2pk B T 



P 2 
2a - 



1 exp 



xp r cj) T (v) exp 



" P 2 

V 

. a 

ma 

k^f^ 



(A13b) 



L/2) 



and <p T (v) is the Maxwell distribution function, Equation 
(J3J). Equations (|Alip and (|A12|) were derived for an ac- 
celerating field towards the right (a > 0). For a < 0, the 
corresponding result for the current is: 



k B T 
2ixm 



pi exp 



-Lip 



k B T 



(A14) 



APPENDIX B 

In section IIVI we constructed the solution of the 
Fokker-Planck equation (|26|l . satisfying the boundary 
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conditions given in Equation (|27|) . The solution consists 
of a linear combination of a non-equilibrium stationary 
homogeneous solution l|28|l . which gives rise to a con- 
stant current, and an inhomogeneous equilibrium state 
(|29|l . which does not contribute to the current. We thus 
found the stationary state (|30|1 , which is of the form 



f(x,v) = Aexp 



max 



knT 



cj> T {v) + 



7 



(Bl) 



where (f> T denotes the Maxwell distribution © . The cor- 
responding current j and density n(x) of the ions are 
given by: 



a 

.1 = ~B 

7 

n{x) — Aexp 



k R T 



D 



(B2a) 
(B2b) 



The imposed incoming ionic fluxes from the thermostats 
determine uniquely the values of A and B (Equations 
(EU and 

When the ion density in the right thermostat is related 
to that in the left thermostat by the Boltzmann factor 



p r = pi exp 



maL 



k B T 



(B3) 



the coefficient B vanishes and the distribution (IB 1|) re- 
duces to the well known equilibrium solution of the FP 
equation 129|) . When relation l|B3|) docs not hold, how- 
ever, equilibrium cannot occur and a stationary current 
will flow through the channel. It can be proved that in 
this situation the ionic number density n(x) is positive, 
but a difficulty arises when one considers the velocity 
distribution in the region of large velocities. In order to 
illustrate the problem let us study the simple limiting 
case of vanishing acceleration. From Equations (|31[) and 
in the limit of a — > we get the asymptotic relations 



A + B = 



Pi 



(B4) 



The evaluation of the point limit (at fixed values of the 
variable v) leads eventually to the distribution 



f a =o(x,v) = 



Pi + Pr (Pi ~ Pr)</' T (0) 
2 1 + 7 L0 T (O) 



(v — 7a;) 4> T (v) 

m 

The distribution l|B7(l is an inhomogeneous solution of 
the the FP equation 

satisfying to the boundary conditions l|27|) . The current 
and the density profile, which is linear, are given by 



J 

n(x) 



Jl -Jr 



1+ 7 L0 T (O) 

Pi + Pr (pi ~ Pr)(t> T (0 ). 

2 1+ 7 L0 T (O) 



(B9a) 
(B9b) 



Whereas n(x) is positive everywhere within the chan- 
nel, the complete distribution f a =o(x,v) is not positive 
definite. For example, when pi < p r , f a =o{x,v) turns 
negative for sufficiently large velocities and thus loses its 
physical meaning. Hence it seems that a physically ac- 
ceptable inhomogeneous stationary state cannot be ob- 
tained from the FP equation. 

Notice that one can define a characteristic velocity jL, 
associated with the finite length channel: ions with veloc- 
ities larger than this will cover the length of the channel 
in a time interval which is shorter than the FP thermal- 
ization time 7 . It turns out that when v < 7L the 
distribution function f a =o{x,v) > 0, and so the difficulty 
appears only for velocities larger than 7L. Clearly, when 
the friction coefficient 7 is large enough the unphysical 
region is quantitatively irrelevant because of the negli- 
gibly small probability weight coming from the Maxwell 
distribution. From a fundamental point of view, how- 
ever, the solution 1|B7(I demonstrates that the FP equa- 
tion is incompatible with the boundary conditions l|27() 
for a channel of finite length. 



(1 + 7 L</> T (0))5 = 7^ T (0) 



Pi 



(B5) 



APPENDIX C 



-(Pi - Pr) 



27ra^ T (0) ' 

The stationary distribution (|Bip takes the form 



lim /(x, v) 



7 



Pi + Pr , T ( \ 

— 2 — 4> ( v ) 

(Pi - Pr) 



(B6) 



27ra0 T (O) 1 + 7L0 T (O) 



x I 4> (y — a /i) — ex P 



k B T 



In this section, we describe the solution of the ki- 
netic equation (|43fl . f° r i° n s flowing through a chan- 
nel under the influence of a uniform accelerating field, 
a Fokker-Planck friction mechanism and stopping traps 
distributed on average uniformly, in the limit where the 
channel is infinitely long and the problem is spatially ho- 
mogeneous. 

The first term on the r.h.s. of Equation (|43f) describes 
the effect of the stopping traps (in terms of a balance 
between loss of particles of dimensionless velocity u and 
gain of those with u = 0). This term may be recast in 
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the form: 



Differentiation of lC6al) leads to 



P^6(u)J + dw\w\F(y,w)-\u\F(y,u)} (CI) 

P d \ f +0 ° 
= —~r- < sgn(u) / dw\u — w\F(w) 

2 du y J_ oc 

r+oa 

— \u\ / dw sgn(« — w)F(w) 



P_d_ 

2du 



sgn(u) / dw w sgn(u — w)F{w) 



u-a+-)F(u) 



where sgn(w) = 9(u) — 9{—u). 

Substituting 1|C1|) into l|43|) . and integrating over u, one 
obtains the following relation: 



(C2) 



= — /3sgn(tt) / dw w sgn(it — w)F(w) + C 

2 J-oo 

The integration constant C on the r.h.s. of this relation is 
determined by the boundary condition \iuy u ^± 00 F(u) = 
0, which shows that 



f +oc 
C = / dw wF(w) 



(C3) 



where j is the dimensionless current. Gathering results, 
Equation l|C2(l may be cast in the form: 



(C4) 



0<0(u) dwwF(w) +9(-u) dwwF(w) 



u-a+—\ F(u) 
du ' 



The structure of the integro-differential equation (|C4ll 
suggests seeking a solution of the form: 

F{u) =6(u)F + {u) + 9(-u)F-{u) (C5) 

where F + and F~ satisfy the equations: 



u-a+ j- ] F + (u) =-(3 j dwwF + (w) (C6a) 



d 2 dF+(v) 
^F+(u) + (u- a)^-^ + (1 - f3u)F+(u) 
du 

= (C8) 



du 2 du 
= 

Seeking a solution of the form 

F+(u) =cxp [-{u-af/A]G + {u) 



(C9) 



and substituting in l|C8() . one arrives at the following dif- 
ferential equation for G + (u): 



d 2 G+(u) 
du 2 



1 (u-a + 2P) 2 
~+P(j3-a) 







G+(u) 
(CIO) 



the solution of which is: 

G+(u) = D m _ a) {u-a + 2!3) (Cll) 

where D p (z) is a parabolic cylinder function. Hence: 

F+(u) = exp [-(« - a) 2 /4]^ (/3 _ Q) ( W - a + 2/3) (C12) 

Proceeding along similar lines, one finds the solution of 
]6bl) in the form: 



F~(u) =exp[-( M - a ) 2 /4] J D (3(/3+Q )(-w + a + 2/3) (C13) 
and hence: 

F(u)= exp[-( U -a) 2 /4] (C14) 
x j^i 6(-u) D w+a) (-u + a + 2(3) 

+A 2 6{u) D m _ a) {u -a + 20)) 



A relation between the coefficients A\ and A 2 follows 
from the boundary condition l|C7(l . namely: 

AiD fj{p+a) (a + 2(3) = A 2 D fj(p _ a) (-a + 2(3) (C15) 

from which the general solution 1441) in the main text 
follows. 



u-a + — ] F~(u) = -P j dwwF-(w) (C6b) 



to be solved, subject to the boundary condition: 

F+(0) -F~(Q) =0 (C7) 
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